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We present a comprehensive study of the thermodynamic properties of the three-dimensional fermionic Hub- 
bard model, with application to cold fermionic atoms subject to an optical lattice and a trapping potential. Our 
study is focused on the temperature range of current experimental interest. We employ two theoretical methods 
- dynamical mean-field theory and high-temperature series - and perform comparative benchmarks to delimit 
their respective range of validity. Special attention is devoted to understand the implications that thermody- 
namic properties of this system have on cooling. Considering the distribution function of local occupancies in 
the inhomogeneous lattice, we show that, under adiabatic evolution, the variation of any observable (e.g., tem- 
perature) can be conveniently disentangled into two distinct contributions. The first contribution is due to the 
redistribution of atoms in the trap during the evolution, while the second one comes from the intrinsic change 
of the observable. Finally, we provide a simplified picture of a recently proposed cooling procedure, based on 
spatial entropy separation, by applying this method to an idealized model. 

PACS numbers: 67.85.-d, 03.75.Ss, 05.30.Fk, 71.10.Fd 



I. INTRODUCTION 

One of the main ongoing efforts in cold atom gases is the 
investigation of strongly correlated phases. Our theoretical 
understanding of strongly correlated phases is, in general, 
far from complete. For example, computational studies, of 
fermions in particular, are severely limited. As a result, the 
possibility of performing analog simulations of model Hamil- 
tonians, using cold atoms in optical lattices^, has raised 
great hopes. The remarkable controllability of cold atom 
systems, allowing, for example, the application of a specific 
time-dependent perturbation, has also opened the possibility 
of studying strongly correlated systems in regimes inaccessi- 
ble to solid state materials, especially away from equilibrium. 

Among all model Hamiltonians relevant to the physics of 
strong correlations, the Hubbard model has attracted the great- 
est attention. On one hand, this model is the simplest to have 
important competition between kinetic and potential energies. 
It also plays in the field of strong correlations a somewhat 
analogous role to the one played by the Ising model in classi- 
cal statistical mechanics. On the other hand, its physics is of 
direct relevance to high-temperature superconductivity, a phe- 
nomenon still veiled in mystery. In fact, this new field often 
dubbed 'condensed matter of light and atoms' was pioneered 
by the theoretical prediction 1 and experimental observation of 
an incompressible regime, characteristic of a Mott insulator, 
and of the transition between this phase and itinerant ones, 
first for bosons 6 and recently also for fermion^H. For recent 
reviews, see, e.g., Refs. |4|5l 

In order to make progress towards the ultimate goal of per- 
forming analog simulations of model Hamiltonians, a good 
synergy between experimental efforts and theoretical investi- 
gations is crucial. For example, theoretical inputs are useful in 
establishing maps, in parameter space, of the location of the 
different phases present in a realistic setup which takes into 



account the trap confining potential. Moreover, experiments 
are currently confronted with the great difficulty of cooling 
fermions in optical lattices to sufficiently low temperatures to 
reach many interesting strongly correlated phases. This rel- 
evant temperature range is significantly lower than the one 
corresponding to mere quantum degeneracy. Theoretical con- 
trol over these issues is greatly needed and requires a quan- 
titative understanding of the thermodynamic properties of the 
Hubbard model, both for the homogeneous system and in the 
presence of a trap. 

In this article, we perform a comprehensive study of 
the thermodynamic properties of the homogeneous three- 
dimensional fermionic Hubbard model, and of cold fermionic 
atoms in a three-dimensional optical lattice subjected to a trap- 
ping potential. We focus on the range of temperature which 
is of direct interest for current experiments as well as for the 
next generation of experiments. Particular emphasis is put on 
aspects related to cooling of the system. 

The main theoretical technique used in the present study 
is dynamical mean-field theory (DMFT). This approach (re- 
viewed, e.g., in Ref. 9) is a controlled approximation which 
is able to capture the competition between the kinetic en- 
ergy, that tends to delocalize atoms over the whole lattice, 
and the repulsive potential energy that prevents atoms from 
occupying the same site, hence promoting localization. We 
also use another theoretical approach, namely high order high- 



temperature series expansions 



One of the goals of the 



present article is to delimit, in parameter space, the respective 
range of validity of each of these approaches, and to provide 
a benchmark for their use through quantitative comparisons. 

This article is organized as follows. In Sec. |IIJ we de- 
fine the model considered in this article, specify notations 
and conventions and briefly outline the theoretical methods 



used. In Sec. Ill we provide detailed results for the thermo- 
dynamics of the homogeneous Hubbard model in three dimen- 



2 



sions. Finally, the effect of the trapping potential is considered 
in Sec. |IV] with several applications geared towards cooling 
fermionic cold atom systems. 



— ln(Z)//3, in (3 J, where Z = Trcxp(— (3H) is the grand 
partition function. The series for the grand potential in a uni- 
form system (V = 0) can be written: 



II. MODEL, METHODS AND CONVENTIONS 



(3fl = In z + J2 (t3J/z ) m X^ (w, C), 



(2) 



Cold fermionic gases with two hyperfine states loaded in 
an optical lattice potential can be described in a wi de range of 
parameters by the fermionic Hubbard mode l 1 ! 14 ! 15 ! . The model 
Hamiltonian is given by 



(1) 



where a is the creation operator of a fermion with spin a at 

site j and hj, a = ct a c - a is the density operator at site j. rj 
is measured in units of the lattice spacing. The spin a =4_, f 
labels the two hyperfine states of the atoms, and labels 
neighboring sites. 

The first term in the Hamiltonian describes the kinetic en- 
ergy of the atoms with J the hopping amplitude between 
nearest-neighbor sites. In this paper, we mainly consider a 
three dimensional cubic lattice and use the half-bandwidth 6 J 
as our energy units. The second term represents the s-wave 
scattering between fermions in different hyperfine states, and 
its strength U is proportional to the s-wave scattering length 
which can be tuned over a wide range using a Feshbach res- 
onance. We often use a grand-canonical description of the 
system in which the filling can be adjusted by the chemical 
potential The last term in the Hamiltonian represents the 
external trapping potential. Aspects related to this trapping 
term will be discussed in more detail later. 

We investigate the properties of the system in a moder- 
ately high temperature regime, typically ksT/QJ > 1/10, 
or /36J < 10 where j3 = l/ksT is the inverse temperature. 
This is the temperature regime of interest in current ultracold 
quantum gases experiments 1 ^. Also one should note that we 
often use fc# = 1. 

At moderate temperature and commensurate filling (one 
atom per site) the homogeneous model displays a crossover 
between a liquid phase at weak interaction and an incompress- 
ible regime at larger interaction. The latter is characterized by 
suppressed density fluctuations and is a Mott insulator. The 
properties of these two phases will be discussed in detail in 
the next sections. A simple picture of these two regimes can 
already be obtained by considering the extreme limits of van- 
ishing interaction and vanishing hopping. In the first case, 
the system only consists of a mixture of non-interacting free 
fermions. The fermions are delocalized and form a Fermi sea. 
In the case of vanishing hopping, the so-called atomic limit, 
the system consists of disconnected sites, and the problem of 
a single site can be solved. The atoms localize and a Mott 
insulator forms at half filling. 

The atomic limit can be seen as the zeroth order term in 
a high-temperature expansion of the grand potential, SI = 



where ( = exp(/3/i) is the fugacity, £1 = Q/N, w — 
cxp(— /3U) and z = 1 + 2C, + C, 2 w is the partition function 
of a single site in the atomic limit. The atom ic limit can be 
improved by systematically calculating x 2i t 13 l 17 H 19 l higher order 
series coefficients, X^ m \ Thermodynamic quantities at high 
temperatures can be accurately derived from leading terms in 
Eq.0 

For intermediate couplings and temperatures, the solution 
of the Hubbard model is highly non trivial. In this article, we 
mainly use dynamical mean-field theory (DMFT) 9 to explore 
this regime. DMFT is an approximate method in which the 
lattice model is replaced by a self-consistent impurity model. 
This model can be solved using highly accurate numerical al- 
gorithms such as the strong-coupling CT-QMC algorithm 20 
which we use in the present work. Although approximate, 
DMFT is a controlled approximation: it is exact in both non- 
interacting and atomic limits, and bridges the gap between the 
two limits. Furthermore, it becomes mathematically exact in 
the formal limit of infinite lattice coordination. Physically, 
DMFT neglects spatial correlations, but treats accurately local 
quantum fluctuations. Besides DMFT , we a lso use high-order 
high-temperature series expansions ' 13 ! 17 ^ . i n f ac t ; one of 
the aims of this article is to provide benchmarking of DMFT 
in the high-temperature domain from comparison with these 
series. More details on the regions of validity of the two meth- 



ods are given in Sec. Ill C 



Finally, at low temperature (/36J > 16), a phase transi- 
tion into an antiferromagnetically order phase occurs^ED, this 
regime will not be considered in the rest of this paper. 



III. HOMOGENEOUS HUBBARD MODEL 

In this section we focus on the properties of the homoge- 
neous Hubbard model. We discuss in particular the density, 
the double occupancy and the entropy in this model. These 
quantities are very important in the characterization of ultra 
cold quantum gases. The effect of the trap will be considered 
inSec.HV] 



A. Density, double occupancy and entropy 

The physics of the competition between kinetic and poten- 
tial energy in the Hubbard model is exemplified by the behav- 



ior of the density per site n. 



l 3,i. 



as a function 



of the chemical potential /.i. For simplicity, since we are con- 
sidering the homogeneous system, we drop the site index in 
this section. At low enough temperature, this quantity has a 
qualitatively different behavior for weak and strong interac- 
tion (Fig. [TJ. At weak coupling (cf. U/6J < 1 in Fig. [TJ, 
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FIG. 1: (Color online) Density n versus chemical potential fj, for 
different interaction strength U and at fixed inverse temperature 
/36J = 10 (DMFT). For large interaction strength a clear Mott 
plateau of filling n = 1 develops. 



FIG. 3: (Color online) Double occupancy d versus chemical potential 
li at fixed inverse temperature /36J = 10 and different values of the 
interaction strength U (DMFT). 




-2 2 4 6 8 

U./6J 



FIG. 2: (Color online) Density n versus chemical potential \i for 
increasing interaction U and fixed inverse temperature /36J = 1 
(DMFT). Only a reminiscent behavior of the Mott plateau is left for 
large interaction strength. 



n(fi) has a smooth evolution from n = for fi <C — 6J to 
n = 2 for fi ^> 6 J + U. In this case, the compressibility 
k = dn/d/j, is always finite when < n < 2. In the opposite 
limit of strong coupling (cf. U/6J > 2 in Fig.[T|), the potential 
energy disfavors the presence of two fermions on a single site 
and n develops a plateau around fi = {7/2. On this plateau, 
characteristic of the Mott insulating state, the compressibil- 
ity k vanishes in the zero temperature limit. As the chemical 
potential must be of the order of U to overcome the potential 
interaction and to generate doubly occupied sites, the central 
plateau widens with increasing interaction. 

Increasing the temperature causes a softening of these fea- 
tures. Fig.|2]shows the density plotted at fcsT/6J = 1. At 
weak interaction strength, thermal excitations reduce the slope 
which characterizes the variation of the density with chemical 
potential. For strong interactions, large temperatures of the 
order of U (the Mott gap) are needed to generate thermal ex- 
citations and to destroy the Mott plateau. One can see com- 
paring Figs.[T|and|2]how the plateau found at ksT/6J = 0.1 
for U/QJ > 2.5 disappears at high temperature only leaving 



behind an inflection in the [7/6 J = 5 curve. 

The double occupancy d = (n>j,inj,t) (Figs. [3] |4jl offers 
another perspective on the Mott phenomenon. This quantity 
is particularly interesting as in cold atom systems it can be 
directly measured^. Hence, understanding how double oc- 
cupancy varies is very useful to characterize the state of an 
experimental system. Furthermore, this quantity is directly re- 
lated to the potential energy of the system which is given for 
a site by U d. In Fig. [3] we display the behavior of d as a func- 
tion of the chemical potential /i. At weak coupling, d tracks 
the behavior of the density (for the non-interacting system, 
d = n 2 /4) and the two quantities contain essentially the same 
information. However, for large coupling, d is strongly sup- 
pressed when the chemical potential is lower than the potential 
energy U, and increases steeply in the region [i > U, precisely 
the region where the density increases above unity. This qual- 
itative difference between weak and strong coupling is well 
apparent on Fig.|4]where we plot d versus n. For U/6J = 0.5, 
the curve is nearly quadratic while for U/6 J = 5 the double 
occupancy d is close to zero for n < 1 and linear for n > 1. 
The behavior of d at strong coupling is to be contrasted with 
the behavior of the density. For strong coupling, n presents 
a plateau at n ~ 1, and hence can be used to determine the 
chemical potentials for which n drops below unity as well as 
for which n exceeds unity. Instead, d is suppressed in the 
whole region /i < U, and is only sensitive to the point at 
which n exceeds unity. This is advantageous when it comes 
to experimentally detecting the Mott stat d 10 l 23 l In the inset of 
Fig. [4] we also display a plot of d/nas a. function of n. 

Another quantity of fundamental importance is the entropy 
as preparation of cold atom systems is often done almost adi- 
abatically. Therefore, the entropy of the system, rather than 
the temperature, is a conserved quantity which can be used as 
a constant to characterize the system evolution. We calculate 
the entropy by integrating a fit to the DMFT energy data start- 
ing from infinite temperature. The high temperature regime 
(/36 J < 1) is approximated using a second order high temper- 
ature expansion. 

In the weakly interacting regime (Fig. [5j, the entropy per 
site s has a simple evolution with temperature and chemical 
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FIG. 4: (Color online) Double occupancy d versus density n at fixed 
inverse temperature /36 J = 10 and different values of the interaction 
strength U (DMFT). Inset: d(n)/n at the inverse temperature /36 J = 
10 and different values of the interaction strength U (DMFT). 




FIG. 5: (Color online) Entropy per site s versus chemical potential fi 
at the weak interaction strength 17/6 J = 0.5 for different values of 
the inverse temperature /3 (DMFT). 



potential: it is maximal at fi = U/2 and decreases mono- 
tonically with increasing \fi — U/2\. With decreasing tem- 
peratures, s(/z) decreases uniformly over the whole chemical 
potential range. 

In the strongly correlated regime (Fig. RS, s(fi) presents 
non-trivial features. Already at relatively high temperatures, 
the entropy has two maxima located approximately at fi = 
and ii = U, while it has a minimum at the particle-hole sym- 
metric point n = U/2. From an atomic limit point of view, 
the presence of these two peaks in the entropy at /i = and 
fi = Uis easily interpreted. These two chemical potential val- 
ues correspond to the regions where the density crosses over 
rapidly between n ~ and n ~ 1 (ji ~ 0), and between n ~ 1 
and n ~ 2 (/i ~ E7). The charge fluctuations in these regions 
are maximal and contribute significantly to the entropy. 

The temperature evolution is also interesting (Fig. [6j: the 
entropy for fi < and fi > U decreases quickly with low- 
ering T, while in the region fi ~ Z7/2 the temperature- 
dependence slows down and the entropy approaches a finite 
value, s(/i = E7/2) — > In 2. This value reflects the fact that, 
in the Mott insulator regime, the system gets frozen into a 



FIG. 6: (Color online) Entropy per site s versus chemical potential 
fj, at the strong interaction strength U /6 J = 5 for different values of 
the inverse temperature j3 (DMFT). The curves legend is the same as 
in Fig. [5] 



local moment state (with two spin-states per site) in a rather 
extended temperature range below the Mott gap. Naturally, 
we expect that, as the system is cooled down further, the en- 
tropy will eventually decrease again below In 2 as magnetic 
correlations develop between local moments. As DMFT ne- 
glects spatial correlations in the paramagnetic phase, within 
this theory, the decrease in entropy only happens right at the 
Neel transition where long-range antiferromagnetic order sets 
in (this Neel transition occurs at a lower temperature than 
the ones studied in the present paper). However, in reality, 
the entropy will start deviating significantly from In 2 abo ve 
the Neel transition, as short-range correlations develorJEEED. 
Nevertheless, the single-site DMFT description can be re- 
garded as a good description of the paramagnetic Mott in- 
sulator down to a characteristic temperature. We will see in 



Sec. Ill C how it is possible to evaluate this temperature with 



the help of high temperature series expansions. 

In Fig. [7] a similar behavior is observed by plotting the en- 
tropy versus the density. For low interaction strength, a max- 
imum of the entropy is found at half filling. In contrast, at 
large interaction strength and intermediate temperatures, a dip 
in the entropy arises at half filling, which is due to the freez- 
ing of the density degree of freedom in the Mott insulator. At 
low temperatures, the entropy decreases linearly in tempera- 
ture for the liquid regimes away from half-filling, while the 
value at half filling remains fixed to In 2 (for a range of tem- 
peratures, see above). 

In Fig. [8] we additionally plot the entropy per particle s/n 
as this quantity will be useful later on to better understand the 
behavior of the entropy in the presence of a trapping poten- 
tial. We note in particular that upon dividing the entropy by 
the particle number, the entropy per particle decreases with 
increasing particle number. 
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FIG. 7: (Color online) Entropy per site s versus density n at differ- 
ent interaction strengths U and inverse temperatures f3 (DMFT). The 
curves legend is the same as in Fig. [5] 
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FIG. 8: (Color online) Entropy per particle s(n) /n at different inter- 
action strengths U and inverse temperatures /3 (DMFT). The curves 
legend is the same as in Fig. [5] 



B. Pomeranchuk effect 

The double occupancy provides a wealth of information on 
the physics of cold atom systems. For example, this quantity 
can be used as a ther mometer in a certain regime of temper- 
atures and interaction g 16 l 22 l 23 l 27 l 28 [ To obtain a good under- 
standing of how this quantity behaves in a trapped system, we 
first show here the dependence of the double occupancy on 
temperature in an homogeneous system (Fig. [9]). 

For moderate interaction strength, the double occupancy 
presents a non-monotonous behavior at low temperatures. The 
appearance of the initial decrease in the double occupancy 
with increasing temperature (from T = 0) is analogous to the 
Pomeranchuk effect observed in liquid Hel ium- 3 a nd has been 
discussed previously in the half filled cas o 15 l 24 l 29 l The reason 
behind this decrease of the double occupancy with increasing 
temperature is that the system prefers, when heated, to local- 
ize the atoms. For this localized state, the (spin) entropy is 
larger than for a state where fermions form a Fermi liquid. In 
this regime, the minimum of d is determined by the quasi- 



particle coherence temperature. Since, in the particle-hole 
symmetric case, the coherence temperature decreases with U, 
the 'Pomeranchuk' effect occurs if £7/6 J is not too large and 
disappears in the Mott insulator. Away from particle-hole 
symmetry, the system behavior is quite different. In that case, 
the 'Pomeranchuk' temperature at which d(T) has a minimum 
increases with U. We also note that the double occupancy has 
recently been studied below and close to the Neel transition 
to the antiferromagnetic ordeP, a regime that we do not con- 
sider here. In the antiferromagnetic phase, the coherent align- 
ment of spins causes an increase of the double occupancy. 

As the entropy and double occupancy are related by the 
thermodynamic (Maxwell) relational 



ds_ 
dU 



dd 



(3) 



the Pomeranchuk minimum in d(T) translates into a non- 
monotonous behavior of the entropy as a function of the in- 
teraction strength. For this quantity, a maximum is found at 
sufficiently low temperatures. As in the case of the double oc- 
cupancy, this effect persists away from half filling (Fig(9]). We 
will see later in Sec. |IV C] the consequences of these consider- 
ations for systems subjected to a trapping potential. 

In the strong interaction limit, the low-energy properties be- 
come similar to the properties of the Heisenberg model. The 
entropy for this model has been recently studied in Ref.l26l 



C. Regimes of validity of DMFT and of high-temperature 
series: a comparative study 

In this section, we perform a comparative study of the valid- 
ity of the two theoretical approximations used in the present 
paper: DMFT and high-temperature series. The motivation 
for doing this is twofold. On the one hand, in regimes where a 
converged result can be reliably extracted from the series ex- 
pansion (which is an exact technique implying no further ap- 
proximation), we can use this comparison to benchmark the 
validity of DMFT as an approximation to the thermodynamic 
properties of the three-dimensional Hubbard model. On the 
other hand, by comparing different orders of the series expan- 
sion, we can delineate the regime in which this method can be 
used reliably and in which regime DMFT is a better option. 
In Fig. 



10 and Fig. 11 



we display the results obtained by 
the two methods for the entropy and the double occupancy as 
a function of temperature, at intermediate interaction strength, 
for two different chemical potentials. Different orders (up to 
order 10) of the series expansions are also compared in these 
figures. It is immediately clear from these two figures that the 
temperature down to which the series expansion can be safely 
used, and consequently down to which a reliable assessment 
of DMFT can be made, depends considerably on the value of 
the chemical potential (or density). For /i/6J = 4 (Fig. 10 1, 
corresponding to a quite high density, the agreement between 
the different orders of the series expansion, and the overall 
agreement with DMFT is essentially perfect over the whole 
range of temperatures considered. In fact, only very small 
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FIG. 9: (Color online) Left column: double occupancy d as a func- 
tion of temperature T for different interactions strength U (DMFT). 
(from top to bottom: n = 1, n = 0.8 and n = 0.6). Right column: 
entropy s as a function of U for different temperatures (DMFT). 
(from top to bottom: n = 1, n = 0.8 and n = 0.6). The inverse 
temperature /3 is measured in units of 1 /6 J. 
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FIG. 10: (Color online) Comparison of the results for tempera- 
ture dependence of the entropy s and double occupancy d at fixed 
chemical potential /i/6J = 4 and intermediate interaction strength 
U/6J = 2.5 obtained by high temperature series and DMFT. 
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deviations of the second order expansion can actually be de- 
tected. This agreement provides strong evidence for the cor- 
rectness of both methods down to quite low temperatures in 
the high-density regime (and by particle hole-symmetry, also 
in the low density regime). In contrast for /x/6 J = 2.45 (cor- 
responding to an intermediate filling n rj 1.3), it is clear from 
Fig. [TT] that the different orders of the series expansion start 
to deviate from one another already at a rather high tempera- 
ture [36 J ~ 4. Below this temperature, the series expansion 
method becomes unreliable. This breakdown happens approx- 
imately at the same temperatures for the entropy and the dou- 
ble occupancy. As shown below, this regime of intermediate 
densities is the hardest one for the series expansion method. 

By comparing the different quantities as a function of the 
chemical potential at fixed temperature, we identify more 



clearly these different regimes of density. In Fig. 12 we show 
the results for the entropy per site s as a function of chemical 
potential at a fixed temperature j36J — 5. Here one can clearly 
identify three different regions. At very low (and very high, 
due to particle-hole symmetry) chemical potential /x, where 
the density of the system is close to zero (close to two parti- 
cles per site), both DMFT and the series expansion give re- 
liable results down to fairly low temperatures (at least of the 
order of (36 J ~ 10). In contrast, in the intermediate region 



FIG. 1 1 : (Color online) Comparison of the results for temperature 
dependence of the entropy s and double occupancy d at fixed chem- 
ical potential /J./6J — 2.45 and intermediate interaction strength 
U/6J = 2.5 obtained by high temperature series and DMFT 
method. 



around /j, ~ and ~ U where the density crosses over from 
to 1 and from 1 to 2 respectively, the series expansion breaks 
down at a rather high temperature. Already at /36J ~ 5 it is 
not reliable anymore. Very close to half filling, the opposite 
trend is found. In this region, the series expansion gives reli- 
able results down to quite low temperatures, actually lower (as 
further detailed below) than the temperature at which DMFT 
ceases to be an accurate approximatiorP} 

In Fig.[T3] we summarize the validity of the high tempera- 
ture series expansion by showing the difference between dif- 
ferent orders. For this plot, it appears that the breakdown tem- 
perature for the series expansion depends sensitively on the 
chemical potential and that the method works best close to 
very low and very high filling, as well as exactly at half fill- 
ing. 



In Fig. 14 we focus in more detail on the comparison of the 
entropy at half filling (n = 1, /i = U/2). As mentioned pre- 
viously, the single site DMFT approximation overestimates 
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FIG. 12: (Color online) Comparison of the results for the entropy 
s versus the chemical potential at intermediate interaction strength 
(7/6 J = 2.5 obtained by high temperature series and DMFT 
method. 
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FIG. 14: (Color online) Entropy per particle in a system at half filling 
and intermediate interaction strength U/6J = 2.5 obtained by series 
expansion, DMFT and QMC (for the Heisenberg modelj^. 
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FIG. 13: (Color online) Contour plot of the absolute value of the 
difference between the results obtained with 10th and 6th order series 
for the entropy versus chemical potential and inverse temperature at 
U/6J = 2.5. 



the entropy at half filling, yielding a value In 2 in the large- 
ly limit because it neglects short-range magnetic correlations. 
This value also corresponds to the mean field value for the 
entropy of the Heisenberg model in its paramagnetic phase. 
In reality, the value of the entropy is reduced by short-range 
antiferromagnetic correlations^. In Fig. 
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this limitation 

of single-site DMFT is exposed as we compare entropy ob- 
tained with DMFT to the ones given by the series expansion 
and the Heisenberg model. While the DMFT curve saturates 
at In 2 at low temperature, the different orders of the series ex- 
pansion reach a lower value of the entropy before diverging. 
The Pade expansion of the high temperature series expansion 
lies slightly above the entropy of the Heisenberg model cal- 
culated in Ref. |26]by QMC simulations. This situation rep- 
resents a case in which the high temperature expansion can 
be used down to a lower temperature than single-site DMFT. 
Finally, we note that f urther extensions of DMFT (especially 
cluster-DMFT methodi^E!) exist which overcome these lim- 
itations of single-site DMFT, and restore the physical effects 
of short-range magnetic correlations. 



IV. TRAPPED SYSTEM 

In present day experimental setups, an external potential is 
usually present to confine the atomic cloud. This confinement 
can be due to different sources, for example to the focusing 
of the lattice laser beams or a magnetic or dipole trap. In 
most cases, a parabolic form is a good approximation for the 
shape of the confining potential at the location of the atom 
cloud. However, many refinements can be considered to ob- 
tain a more precise spatial dependence of the confining poten- 
tial V(r) (last term of Eq.[TJ. 



A. Local density approximation and local occupancy 
histogram 

Within the local density approximation (LDA), the prop- 
erties of the trapped system at a certain position rj are as- 
sumed to be those of the homogeneous system with the chem- 
ical potential set to the value /tt(fj). LDA has been found to 
be a good approximation for local quantities in a three dimen- 
sional fermionic gas in an optical lattice®^, and was also val- 
idated in the high temperature regime applicable to ongoing 
experiment ;] 10 ! 16 ! However, as LDA neglects the influence of 
surrounding sites with different densities, the main inaccura- 
cies occur close to phase boundaries, where proximity effects 
occur, and also whe n on e computes non-local longer range 
physical observablesSEl 

Using LDA and taking the continuum limit, it is possible 
to describe the system using rescaled variables that do not de- 
pend explicitly on the strength of the confining potential. To 
illustrate this we consider a radial symmetrical form for the 
trapping potential V(r) = Vt(\r\/a) a in d dimensions (a is 
the lattice spacing). Thus, the position-dependent chemical 
potential becomes 



(4) 



with rj the d-dimensional position vector labeling each lat- 
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tice site and /xq the chemical potential at the center of the 
trap. Within LDA, any local observable 0(r\/) in the inho- 
mogeneous system is related to its homogeneous counterpart 
by 0(rj) = Oh(n(rj)), where Oh denotes the corresponding 
quantity in the homogeneous system (in the following we will 
drop the label 'h'). 

The number of atoms, which corresponds to the sum of the 
local occupancy over the whole system, can be expressed in 
the following way: 



N 



d-l 



drr n{r) 



(5) 



with Od-i the surface of a sphere in d dimensions. Changing 
variables to an integration over the chemical potential using 
(HI finally leads to: 



' V t \ d/a 

" - v is) 



d-l 



dp (jj, - fi) °> 1 n(p). (6) 



In this expression ~p = p/6J is the dimensionless chemi- 
cal potential (similarly JL ). The LDA approximation enters 
in the assumption that n is a function of r only through its 
dependence on p. This formula can be easily generalized 
to the case for which the strength of the confining potential 
is different along the different Cartesian directions (which 
is usually the case in experiments) by replacing V t with a 
proper averaged quantity Vt- For example, in three dimen- 
sions Vt = {Vt tX Vt,yVt.z) 1 ^ 3 - Eq. |6] shows that the dimen- 
sionless combination p = iV(Vt/6J) d/a is a quantity that 
does not depend on the strength of the confining potential, 
and hence can be used to describe properties of experimen- 
tal systems regardless of the particular realization of the trap. 
For the one-dimensional case, see Ref. [36] It also shows that 
in this approximation the key quantity is the observable in the 
homogeneous system, e.g., and that everything can be 

derived from it. Obviously the same holds true for all the other 
quantities that can be expressed as a sum of a local quantity 
over the whole system. 

Another very useful way to express averages of observables 
over the trap is to introduce the distribution function of site 
occupancies in the system, defined for < n < 2 as: 



(7) 



and the related quantity in which each site is weighted by its 
occupancy: 



Q(n) — 7ij5(rij — n) = nP(n). 



(8) 



Using again a continuous-space integration over the system 
and changing variables in favor of the chemical potential, one 
obtains: 



q(n) ee (V t /6jf/°Q(n) = ^ " _ p^f^ 1 



a K,(n) 



In this expression, k(ti) is the (dimensionless) compressibility 
of the homogeneous system: 



n{n) 



dn 
dp 



(10) 



fi=fi(n) 



Hence, the local chemical potential and local compressibility 
entirely determine the distribution of local site occupancies in 
the trap. The distribution Q(n) and rescaled distribution q(n) 
obey the sum-rules: 



dnQ(n) = N 



dn q{n) = p 



(11) 



with N the total atom number and p = (V t /QJ) d/a N its 
rescaled form. 

These distribution functions allow to express the average of 
any observable over the trap, within the LDA approximation, 



o(j) = / dn P(n)o(n) = / dnQ{n) 



o{n) 



O EE ^ 

(12) 

where o is the local operator corresponding to the observable. 
Because Q{n) is a normalized distribution obeying the sum- 
rule ( [TT] ), the last expression is particularly useful. As we shall 
see later, when varying an external parameter, it allows to sep- 
arate in a simple manner the changes in O which are due to a 
redistribution of the particles in the trap (reshaping of Q(n)) 
from the contribution due to the intrinsic dependence of the 
local observable on the parameter, already present in the ho- 
mogeneous system. 

In the following, we make use of this description, in com- 
bination with thermodynamic relations, in order to discuss 
cooling or heating of the trapped system as the coupling 
is changed. We will concentrate on the case of a three- 
dimensional lattice in a harmonic potential i.e., d = 3, a = 2. 



B. State diagram 

One consequence of the presence of an inhomogeneous 
trapping potential is that different quantum phases can spa- 
tially coexist in the gas. This can actually be seen as a fa- 
vorable situation, in which several different physical regimes 
can be studied in a single experiment. In an optical lattice 
realizing the three-dimensional Hubbard model, coexistence 
between liquid and Mott insulating regions in the trap were 
for example documented in theoretical studie d 23 ! 34 ! 



In Fig. 15 we display the different regimes expected in 
a three-dimensional optical lattice confined into a parabolic 
trap, as a function of the coupling U /6 J and of the scaled 
particle number p. Different temperatures in the currently ac- 
cessible range are considered. At still lower temperature (not 
displayed), antiferromagnetic long-range order 37 will occur in 
the regimes with a commensurate Mott plateau. The state di- 



(9) agram of Fig. 15 which generalizes to different temperatures 
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the results of Ref.|23] was obtained on the basis of the theoret- 
ical calculations described in the previous section for the ho- 
mogeneous model, using LDA approximation. In addition, in 
Appendix [A] we introduce a simple approximation which al- 
lows to obtain analytical expressions for the various crossover 
lines of the state diagram in the low-temperature regime. 

The state diagram displays four characteristic regimes (la- 
beled L, B, Mc and Ms), which are illustrated by the four cor- 
responding density profiles n(r) and local occupancy distri- 
bution functions q(n) calculated at four representative points 
and displayed in Fig. [16] 



For low interaction strength (regime 'L', Fig. 16 a) the den- 
sity profile adjusts to the trapping profile and the system re- 
mains a Fermi liquid everywhere in the trap. With increasing 
temperature the density distribution broadens. The weighted 
particle number distribution q(n) displays a maximum at fill- 
ing unity which, according to Eq. Q reflects the smaller com- 
pressibility at that filling. A rather sharp drop is seen towards 
larger occupancies which represents the center of the trap, 
whereas a slower decay occurs towards smaller particle num- 
bers, due to the tails of the density distribution. Increasing 
the temperature shifts weight from larger occupancies towards 
smaller occupancies. 

For very large values of the scaled particle number p, a band 
insulator with n = 2 forms in the center of the trap (regime 
'B', Fig. 16 c). The pinning to n = 2, and hence the band 
, is destroyed by increasing the temperature. In the 



Fig. 
insulator 



presence of a large band insulating region, the corresponding 
distribution q(n) displays a sharp peak at filling n = 2. In- 
creasing the temperature, this peak decreases and the weight 
moves to lower occupancies. 

For larger interaction strength (regime 'Mc', Fig. 16 e) a 



Mott-insulating region appears, in which the density is pinned 
to n = 1 particle per site. Close to the boundary of the trap, 
the Mott insulating region is surrounded by a liquid region. 
The Mott-insulating region shows up in q(n) as a large and 
narrow peak at filling n = 1 with a sharp edge on the large 
occupation numbers side. The peak reflects the essentially 
vanishing compressibility of the Mott insulator. Increasing the 
temperature decreases the size of the Mott insulating plateau, 
and results in a shift of the weight from filling one to low 
densities. 

Increasing the number of atoms in the trap at large interac- 
tion strength can increase the pressure exerted on the atoms, 
and can cause the occurrence of a liquid region with filling 
larger than one in the center, surrounded by a shell of Mott in- 
sulator with n = 1 (regime 'Ms', Fig. 16 g). Correspondingly, 



the sharp peak in q(n) broadens somewhat. 

Recently, experimental evidence of the Mott insulating re- 
gion has been reported^. This has been achieved by ob- 
serving the suppression of the double occupancy in the Mott- 
insulating regiorPDSm an( j the compression of the cloud as a 
response to the variation of the external trapping potential 8 . 

We note that in bosonic two-dimensional gases the density 
profiles n(r) and therefore the occupancy number distribu- 
tions q(n) can nowadays be measured with a very high spatial 
resolutio n 38 " 40 ! In three dimensional gases, the integrated col- 
umn density can be measured for example by using an elec- 
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FIG. 15: (Color online) State diagram of the gas in a three- 
dimensional optical lattice with parabolic trapping, for different tem- 
peratures (DMFT). The four characteristic regimes (see text) are la- 
beled by: B (band insulator in the center of the trap), Mc (Mott in- 
sulator in the center of the trap, shaded areas), Ms (shell of Mott 
insulator away from the center) and L (liquid state). For each temper- 
ature the (crossover) lines indicate, from bottom to top, the p values 
at which the central density takes the values 0.995, 1.005 and 1.995. 
The gray dashed line marks the crossover from the liquid to the Mott 
state with increasing interaction. The crosses indicate the points for 
which the density profiles are plotted in Fig.|16| 



tron microscope 41 . Furthermore, new techniques such as the 
immersion of a single trapped ion into the atomic gas are being 
developed to locally measure the density in three dimensional 
systems as well 4 ^. 



Temperature changes in the trap during an adiabatic 
evolution 



Cold atom clouds are almost perfectly isolated from their 
environment. Therefore, assuming that manipulations can be 
performed adiabatically, the quantity that is conserved during 
the evolution of the system is the entropy. However, the tem- 
perature will in general change. Consequently, studying the 
effects of an isentropic change of parameters on the system 
temperature is important. Here we focus our attention on the 
increase of the interaction strength, U, and identify if such a 
change can help to reach lower temperatures in a trapped sys- 
tem, i.e. whether cooling occurs. 

The change in temperature induced by a change in a certain 
parameter x at constant entropy is ST/5x\s- If this quantity is 
negative (positive), an increase of the parameter is associated 
with cooling (heating). As we now show, one can consider 
alternatively the change in entropy at constant temperature 
8S/8x\t- Indeed, the location T(x) of isentropic lines in the 
(x, T) plane is defined by the equation S[T(x); x] — const.. 
Taking a derivative of this equation, we obtain: 



ST 

Sx 



55 
ST 



ss 

Sx 



0. 



(13) 



Denoting by c = TSS/ ST the specific heat of the system, we 
finally obtain the expression of the relative cooling rate in the 
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FIG. 16: (Color online) Density profiles (left column) and occupancy 
distribution q(n) (right column) for four typical points in the state 
diagram (obtained with DMFT): a and b (in regime 'L'): U /QJ = 1, 
p = 10. c and d (in 'B' for low-T): U/QJ = 1, p = 20. e and f (in 
'Mc'): U/QJ = 3, p = 10. g and h (in 'Ms'): U/QJ = 3,p= 20. 
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FIG. 17: Schematic evolution of the temperature as a function of a 
parameter x along isentropic lines, in a case where cooling occurs as 
x is increased. The continuous line corresponds to a lower entropy 
per particle than the dashed line, s\ < S2. The figure illustrates that 
for this situation, increasing x while keeping T constant will increase 
the entropy (Eq.|14|l. 



0.25 




FIG. 18: (Color online) Temperature dependence of the number of 
doubly occupied sites D(T) for the interaction strength U/QJ = 1 
and at different scaled particle numbers p = 5, 10, 15 (DMFT). 
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Since c is a positive quantity, we see that there is cooling 
(heating) when ff | T is positive (negative). This is illustrated 
schematically in Fig. [XT] 

In the following, we consider the temperature change un- 
der an adiabatic increase of the coupling U. We observe that 
the derivative |^ | T is related to a derivative of the total dou- 



ble occupancy D = . (n^ p 

relational 



through the Maxwell 



SS_ 
5U 



N,T 



6D 



(15) 



N,U 



So that the relative cooling rate reads: 



1 ST 
f SU 



S,N 



1 SD 

c Jr 



(16) 



Hence, when the derivative of D with respect to temperature 
is negative (positive) there will be cooling (heating) upon an 
isentropic increase of the interaction strength. One advantage 
to use the change of D is that this quantity can be measured 
quite accurately in present experiment*^! 

In Fig. T8p0 we plot D as a function of temperature 
for different interaction strengths and particle numbers. For 
U/6J = 1 (Fig. 18 1, D(T) is a decreasing function of tem- 
perature for all particle numbers. This implies that in the 
weak-coupling regime an increase of U generates cooling. In- 
creasing the interaction the situation gradually changes. For 
[7/6 J = 2 (Fig. [19]), D(T) becomes much flatter and cool- 
ing is restricted to large particle numbers or low temperatures. 
Finally, for U/QJ = 3 (Fig. 20 1, the tendency inverts and at 



N,U 



high temperature, heating occurs. It has to be stressed that the 
absolute value of the derivative drops with increasing interac- 
tion and that at large U the heating or cooling is essentially 
negligible. 

In order to better understand the origin of the cooling or 
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FIG. 19: (Color online) Temperature dependence of the number of 
doubly occupied sites D(T) for the interaction strength (7/6 J = 2 
and at different scaled particle numbers p = 5, 10, 15 (DMFT). 




FIG. 20: (Color online) Temperature dependence of the number of 
doubly occupied sites D(T) for the interaction strength (7/6 J = 3 
and at different scaled particle numbers p = 5, 10, 15 (DMFT). 



heating, we write the derivative 5D/5T\n .u using Eq. (|9jl 
under the form: 



5D 



N,U 



dn 



[q(n)d T d\ n +dd T q(n)} 



(17) 



The usefulness of this expression resides in the clear separa- 
tion of two different contributions. The first term in the in- 
tegrand takes into account changes to the total double occu- 
pancy in the trap due to the intrinsic temperature dependence 
of the double occupancy d(n, T, U) in the homogeneous sys- 
tem. Cooling can occur from this term whenever the Pomer- 
anchuk mechanism discussed above applies. The second term 
instead represents the contribution due to the redistribution of 
the atoms in the trap upon a temperature change. The two 
terms can be calculated separately in order to determine the 
most relevant mechanism behind the cooling observed for low 
interaction. 

In Fig. 21 we plot the first ("intrinsic", dotted lines) and 
second ("redistribution", continuous lines) terms of Eq. 17 at 
(7/6 J = 1 for different atom numbers p. Note that, for read- 
ability, these quantities are plotted versus inverse temperature 
f3 = I/UbT. We notice that the "redistribution" term is al- 
ways negative, and hence the reshaping of the density pro- 



file always induces a cooling effect. In contrast, the "intrin- 
sic" term is positive at high temperature and becomes negative 
only for /36J > 5 (k B T/6J < 0.2). Hence, we conclude that 
at high temperature the cooling is dominated by the redistri- 
bution of atoms in the trap, while at lower temperatures, both 
the redistribution and the intrinsic 'Pomeranchuk' effect con- 
tribute on comparable footing. The latter may even dominate 
at still lower temperatures (e.g., ksT/QJ < 1/8 in Fig. 21 1. 

A qualitative understanding of the behavior of each term in 
Eq. ( fTT) can be achieved from the inspection of the properties 
of d(ri) I n and q(n). The first observation is that d(n)/n is a 
monotonically increasing function of n (cf. inset of Fig. |4]). 
Secondly, from Fig. 16 3 and d we notice that drq(n) (possi- 
bly with the exception of the region around n = 1) is always 
negative for n larger than a certain value n and always posi- 
tive for n smaller than n. Furthermore J dn q(n) = p implies 
that J dn drq(n) = 0. Combining these two observations, we 
conclude that the second term in Eq.([T7]i is generally negative. 
Hence, the redistribution of the atoms indeed produces cool- 
ing in general, as observed above. The presence of the peak 
in q(n) around n = 1 might undermine the reasoning but this 



is never the case for the parameters considered here (Fig. 21 



The behavior of the first term in Eq. 17 is closely connected 
to the Pomeranchuk effect in homogeneous systems. As we 
saw in Fig. [9] in the homogeneous system the Pomeranchuk 
effect is active only at sufficiently low temperature. Therefore, 
the intrinsic contribution in the trap can only lead to cooling at 
low temperatures. This is in good agreement with the results 
inFig.[2T] 

The situation is only quantitatively different at stronger 
coupling. For (7/6 J = 3 (Fig. 22 1, the redistribution of atoms 



still corresponds to cooling, although the absolute value of the 
contribution is roughly an order of magnitude smaller than 
for weaker interaction strength. However, in this case, the 
'intrinsic' contribution becomes dominant and its sign corre- 
sponds to heating (opposite to the Pomeranchuk effect). At 
lower temperature, the intrinsic term becomes negative again 
(cooling) but anyhow the cooling rate in this regime is fairly 
small. 

The conclusion of this section is that in trapped systems an 
increase of the interaction is accompanied by cooling for in- 
teractions weaker than the interaction needed to have a Mott 
insulator. On the other hand, at larger interactions, there is no 
substantial cooling associated with an increase of (7, and even 
slight heating can occur at high temperature. On the whole, 
the dominant contribution to cooling is usually the redistribu- 
tion of the atoms in the trap, although the Pomeranchuk effect 
(intrinsic contribution) can become operative at low tempera- 
ture. 



D. Cooling by trap shaping 

Reaching sufficiently low temperatures to observe complex 
quantum phases is one of the main challenges currently faced 
by cold atom physics experimentalists. In this section, we 
show that by adiabatically reshaping the confining trap, to di- 
vide the system into entropy rich and poor regions, a gas can 
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FIG. 21: (Color online) The two contributions to SD/ST for 
U/6J = 1 (DMFT). Continuous line (left side, top to bottom: 
p = 5, 10, 15, 20): redistribution of atoms in the trap (reshaping, 
second term in Eq. jTT}). Dotted line (right side, top to bottom: 
p = 20, 15, 10, 5): intrinsic change of the double occupancy with 
temperature (first term in Eq.(| 1 7|>). 
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FIG. 22: (Color online) Different contributions of the intrinsic and 
the reshaping term to SD/ST for U/6J = 3 (DMFT). Continuous 
line (left side, top to bottom: p — 5, 10, 15, 20) is the reshaping, 
dotted (left side, top to bottom: p — 20, 15, 10, 5) is the intrinsic 
double occupancy change with temperature. 
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FIG. 23: (Color online) Occupation number (dashed line, left upper 
panel), entropy per particle (solid line, left upper panel), and poten- 
tial profile (solid line, lower panel) for the idealized trapping poten- 
tial. Q(n) (solid line (full system), dashed line (core region), right 
upper panel) as a function of the density. The core region is taken to 
be cylindrical with rn = 51a, the storage region is a rectangular box 
of size 300 x 300 x 100a 3 . Each region has a homogeneous density 
Hcore = 2, n slorag e = 0.31. The average entropy per particle for the 
total system is |^ = 1.65 and = for the core region. The 
ratio of particles in the core region versus the total particle number is 
^ = 0.39. 



into a reservoir but on creating spatially distinct regions of 
high and low entropy that can be isolated from one another, 
and on subsequently removing the high entropy region. Once 
this procedure is completed the remaining system has a much 
lower average entropy per particle allowing for the study of 
interesting phenomena requiring lower temperatures than pre- 
viously attainable. To demonstrate our cooling procedure, we 
use a twofold approach. We first present our method using an 
idealized setup, and in a second time, we revisit with a new 
perspective the experimental setup presented in Ref. |43] We 
further highlight the differences between our cooling method 
and another one recently proposed in Ref. [5T| 



be cooled down by one order of magnitude lower than cur- 
rently achievable using state-of-the-art techniques 43 . Before 
presenting our method, we would like to point out that em- 
ploying adiabatic changes to cool down gases was proposed 
in related contexts. For example, a Bose-Einstein conden- 
sate was experimentally produced by adiabatically deform- 
ing the external trapping potential to increase the gas phase 
space densitjEHIl Furthermore, reshaping the underlying trap 
to create entropy rich regions that are later isolated from the 
remaining sy stem was proposed for bosons loaded into an 
optical lattic cr 6 l 47 [ For fermions confined to an optical lat- 
tice, cooling could be achieved by immersion into a bosonic 
batr l 48 ! 49 ! In Ref. |49j it was cleverly proposed to reduce the 
entropy of lattice fermions in contact with a bosonic reservoir 
by compressing them into a band insulator or more generally 
a gapped phase. This last cooling method requires a trans- 
fer of entropy between two distinguishable quantum gases, a 
process experimentally demonstrated in Ref.l50l 

Our cooling method does not rely on immersing the atoms 



1. Idealized realization 

As explained above, we first present our cooling method us- 
ing an idealized setup uncluttered with experimental details. 
In the following, we show how, given an idealized two-fluid 
model, one can start with a gas, loaded in an optical lattice and 
confined to an harmonic trap, having an entropy per particle, 
e.g., around | initial = 1.7 and obtain a final system character- 
ized by a near zero entropy per particle while keeping about 
40% of the atoms. 

As we have seen in Fig. [8] at high temperatures, the en- 
tropy per particle is largest for low densities. Therefore 
to segregate the entropy in our system, we would like to cre- 



ate two distinct regions (Fig. 23 1: (i) a "core region" with a 



deep trapping potential in which the density nc is close to 
two particles per site, i.e. nc ~ 2 and (ii) a "storage region" 
with a very flat trapping potential in which the density n$ is 
very low, i.e. rig <C 1. In the language of the Q(n) distribu- 
tion introduced in the previous sections, this idealized setup 
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results in two sharp peaks: one at a very low density due to 
the storage region and another one very close to n = 2 due to 
the core. 

To obtain these two regions of unequal densities, we start 
from a gas confined to the usual harmonic trapping potential 
with total entropy St, and adiabatically deform the trap to 



reach the trapping profile presented in Fig. 23 While the total 
entropy remains constant under this deformation, the entropy 
is now inhomogeneously distributed. The entropy per particle 
in the core region will be very low, ideally zero, whereas in 
the storage region the entropy per particle will be quite large. 
In contrast, the temperature Tq remains equal throughout the 
system and is set by the constraint requiring entropy conser- 
vation, i.e. 



S = dn< Qc(n) 



s(n,T ) 



+ Qs(n)^^\ = S T 
n J 

(18) 

where Qc{n) = ^iec 7 ^ ^( n * ~ n ) and Qs(n) = 
Si es n i 3( n i ~ n )- F° r the idealized two-fluid model pre- 
sented here, Eq. [18] simplifies to 

S = Nc s J^M + Ns ±^l = St (19) 



n c 



where Nq and Ns are the number of particles in the core and 
storage regions respectively. 

In this idealized situation an infinitely narrow barrier is used 
to separate the two regions adiabatically. After the separation, 
the high entropy region can be removed and the core region 
can be used to perform the experiment. As the density is uni- 
form through out the core region, the average entropy per par- 
ticle at the time of separation is given by = s ^ ra ^ T °^ . 
Therefore, the final average entropy per particle in the core re- 
gion, ^ fflL , can be much lower than the initial average entropy 
per particle. It is important to note that after this separation 
further adiabatic changes of the system parameters can lead to 
temperature changes, but the average entropy per particle will 
remain constant. 

The minimal entropy which can be reached depends on dif- 
ferent quantities as seen from Eq. 19 As expected the initial 
entropy St is an important factor. Furthermore, the cooling 
procedure is more efficient if the entropy per particle in the 
storage region is much larger than the entropy per particle in 

s(nc,T ) „ s(n s ,T 



the core region, i.e. 



For fixed temper- 



nc ns 

atures, a typical behavior of the entropy per particle with the 
density is shown in Fig. [8] As one can see the decrease in en- 
tropy per particle with increasing density is more pronounced 
at larger temperatures. Hence the procedure will become less 
efficient at lower initial temperature. 

We also find that the differences in entropy per particle are 
largest between very low and very high densities, such that 
the procedure would be most efficient if such densities were 
used for the storage and core regions respectively. However, 
it is important to note that it is not essential that the state sta- 
bilized in the core region is gapped as the band insulator is. 
What really matters is that a sizable difference is achieved be- 
tween the entropies per particle characterizing the densities 
of the core and storage regions. By tuning the trap shape the 
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FIG. 24: (Color online) Occupation number (dashed line, left upper 
panel), entropy per particle (solid line, left upper panel), and poten- 
tial profile (solid line, lower panel) in the presence of the dimple and 
barriers, as a function of the transverse coordinate. The potential is 
offset such that V = is at the bottom of the dimple. Q(n) (solid line 
(full system), dashed line (core region), right upper panel) as a func- 
tion of the density. We chose the following experimentally realistic 



parameters: 



0.5, 



1.8104, Y 



50, 



6, Tb = 15a, 



5 a, 



= 15, u)d = 15a, and 12 • 10 4 atoms. The average 
entropy per particle in the total system is Mf- = 1.95 and in the core 



= 0.198. The ratio of particles in the core region versus 



region 

the total particle number is 
also Fig. 2 of Ref. [43] 



0.404. Obtained with DMFT See 



number of atoms in the core and storage regions can be ad- 
justed. Ideally one would like to create a very large storage 
region with a lot of atoms at very low density, since there the 
entropy per atom would be maximal. However, this situation 
can only be achieved within a certain range due to experimen- 
tal limitations. In particular the trap can only be engineered 
within a certain spatial extension and the final number of par- 
ticle (iVfinai = Nc) should be reasonably large in order to 
perform the actual experiment afterward. 

Before turning our attention to the experimentally realiz- 
able setup, one important comment is in order. We want to 
emphasize that the adiabatic isolation of the core region from 
the storage region is favorable or even necessary to achieve 
cooling. Without proper isolation, the entropy of the storage 
region can flow back into the experimentally relevant region. 
At best this backflow of entropy may simply render the cool- 
ing scheme less efficient, but in the worst scenario it may heat 
up the region of interest. At low temperature, heating may 
occur because the experimentally relevant phase may accom- 
modate more easily the entropy than its parent high density 
state (in our case the band insulator). Therefore, while chang- 
ing back the trap shape to generate the experimentally relevant 
phase, the entropy may flow back into the region of interest if 
the core and storage regions are not properly isolated from 
each other. 



2. Realistic realization 

Let us now turn our attention to the experimentally realiz- 
able trapping potential. To be truly relevant, this trapping pro- 
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file should be achievable using present technology and should 
not require fine-tuning of the system parameters. Such an ex- 
perimentally realistic trapping potential has been described in 
Ref. [43] Here we only briefly summarize the setup. In or- 
der to achieve the required entropy modulation, the poten- 
tial has to allow for a tight trapping in the core region, for 
a wide shallow ring in the storage region and for high poten- 
tial barriers isolating these two regions from each other. To 
produce this profile, we envision to use three elements, (i) a 
shallow harmonic trap (either magnetic or optical), (ii) a dim- 
ple which confines atoms to a small region around the trap 
symmetry axis and helps to create the band insulator, and (iii) 
a cylindrically-symmetric potential barrier to isolate high and 
low entropy regions. The dimple (ii) and potential barrier (iii) 
are produced by red- and blue-detuned laser beams respec- 
tively, creating attractive or repulsive dipole potentials. The 
dimple has a Gaussian profile, while the barrier needs to be 
a narrow annulus. Experimentally this can be realized either 
by setting a tightly focused laser beam in rapid rotation, or 
by engineering the beam profile using phase plates or other 
diffractive optics 52 . Consequently, in addition to the lattice 
potential, the trapping profile is given by 

armomc + Vdimple + H 

arner 

With Harmonic (r) = V h {x 2 + y 2 + "( 2 Z 2 ) / 2 , 

^dimpie(r) = -V d exp (~2(x 2 + y 2 )/wj). 
Karrier(r) = V b exp {-2(y/x 2 + y 2 - r b ) 2 /w 2 ), 

where V{h,d,b} are the potential amplitudes, 7 is a measure of 
the anisotropy of the harmonic trap, u>{d,b} are the waists of 
the Gaussian laser beams forming the dimple and barrier, r& 
is the radius of the cylindrical barrier, and a the lattice spac- 
ing. An example of the resulting trapping profile is shown in 



Fig. 24 The gain in entropy per particle for this particular case 
is approximately one order of magnitude while keeping about 
50% of the atoms. In Fig. 24 in addition to the density distri- 
bution and entropy per particle for the corresponding trapping 
profile, one can see on the right panel from the Q(n) distribu- 
tion that the two sharp peaks of the idealized setup have been 
replaced with broader features. In particular, the lower peak 
is not well defined and removing the storage region results in 
an almost complete extinction of its weight. 

However, experimentally removing the storage region is not 
straightforward due to the presence of the optical lattice po- 
tential. In Ref. [43] several possible methods have been pro- 
posed. Fortunately, since the publication of this article, dif- 
ferent techniques aiming at locally addressing the atoms have 
been developed and could be employed to remove the storage 
region. One of these methods makes use of an electron micro- 
scope and could be used to blast away the storage region in 
a very controlled manned. Another possibility would be to 

use locally trapped ions to remove the atoms from the storage 

• J42I 
region . 

All previous considerations assumed the system evolution 
to be perfectly adiabatic. In a real experiment this will not be 
the case as the deformation of the trap needs to be performed 
within a finite time. In Ref. |43j we showed that, for a one- 
dimensional system, reshaping the trapping potential can be 



performed within an experimentally realistic timescale while 
inducing heat on a scale more than ten times smaller than the 
antiferromagnetic exchange coupling. 

Consequently, using the scheme presented above, exper- 
imentalists could cool down quantum gases to one order 
of magnitude lower than presently achievable while keep- 
ing about half the atoms in the system. Cooling into highly 
sought-after quantum phases could thus be achieved. 



V. CONCLUSION 

In this article, we performed a detailed study of the thermo- 
dynamics of the three-dimensional fermionic Hubbard model, 
for a rather wide range of couplings. We mainly focused on 
the temperature regime fceT/6 J > 0.1, of current interest for 
experiments on cold fermionic atoms in optical lattices. 

Our theoretical study is mainly based on single-site dynam- 
ical mean-field theory (DMFT), a well established theoret- 
ical method based on a controlled approximation in which 
non-local correlations are neglected but local quantum fluc- 
tuations are treated accurately. In addition, we used high- 
temperature series expansion. A comparative study between 
these two methods allowed for a precise assessment of their 
respective range of validity. DMFT is found to be accurate 
down to fairly low temperatures when not too close to half- 
filling (one particle per site). Because of their convergence 
properties, the series expansion are most useful at low density 
or exactly at half-filling. At or close to half-filling, single- 
site DMFT can be trusted only down to T ~ J. Below this 
temperature, short-range magnetic correlations set in, which 
require the use of cluster extensions of DMFT. Hence, our 
study validates the use of single-site DMFT for understanding 
experimental results on cold fermionic atoms in a 3D lattice 
in the currently accessible temperature range, while future ex- 
periments at lower temperature will require extensions of the 
method. 

We considered the implications of the thermodynamic 
properties of the homogeneous system for fermionic atoms in 
an optical lattice confined in a parabolic potential, within the 
local density approximation. A state diagram was established, 
with different regimes for the density profile in the trap. Spe- 
cial emphasis was devoted to the distribution function of site 
occupancies in the trap, a quantity which can be experimen- 
tally measured by imaging techniques with single lattice-site 
resolution. Such techniques have recently become available 
for two-dimensional systems. 

This distribution function proved to be of particular use 
when discussing how a given observable (e.g., temperature) 
changes in the trap as the system evolves in an adiabatic man- 
ner when a parameter is varied. Indeed, it allows for a clear 
separation of two contributions, one corresponding to the re- 
distribution of atoms in the trap during the evolution, and an- 
other from the intrinsic change of the observable. We applied 
these considerations to the temperature change in the trap un- 
der an adiabatic variation of the coupling, and identified the 
regimes where a significant cooling can be expected. How- 
ever, we would like to point out that for these results to be 
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completely applicable to trapped cold atoms confined to opti- 
cal lattices, one would need to take into account several pro- 
cesses that could prevent a fully adiabatic evolution. Consid- 
ering the effects of processes such as lattice heating will be 
the subject of future works. 

Finally, we elaborated on a previously proposed cooling 
mechanism, based on the separation of regions of small and 
large entropy. This procedure is promising for cooling further 
cold fermionic systems by approximately one order of magni- 
tude, a necessary step and major future challenge for access- 
ing and studying experimentally strongly correlated phases in 
those systems. 

Note added: When the present work was near completion 
we became aware that another study, conducted by Fuchs et 
ai, on the thermodynamics of the 3D Hubbard model was 
also in its final stage. The focus of Ref. [53] is on the char- 
acterization of the low temperature region near the antiferro- 
magnetic transition. In the temperature range where our two 
works overlap, the results are in agreement. 
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Appendix A: A simple analytical approximation 

In this appendix, we present a very simple approximation 
that allows us to obtain the state diagram (Fig.[l5]and Ref.l23t 
of trapped fermionic atoms in an optical lattice with very lit- 
tle computational effort. It also allows us to obtain approxi- 
mate analytical expressions for the crossover lines in this state 
diagram, hence providing qualitative understanding into the 
numerical results presented above in this article. 

This approximation is based on the following approximate 
form of the one-particle spectral function of the homogeneous 
system: 



A(w) 



U-A+2D 

1 

U-A+2D 





for uj + /i < 
for - D < uj + 



D 



u < - ^ 

f* 2 2 
+ (1 < f + f 



for 



<lj + [i<U + D 



forw + n > U + D. 



(Al) 

In this expression, D = 6 J is the half-bandwidth and A is 
a parameter that plays the role of the Mott gap (see below). 
The rationale behind this expression is the following. When 



A = 0, (weakly interacting regime) it describes a density of 
states broadened by interaction effects. When A 7^ 0, it de- 



scribes two Hubbard bands (Fig. 25 1 separated by a Mott gap 



The contribution of quasiparticle states appearing in between 
the two Hubbard bands are neglected, because the tempera- 
tures considered in this article are typically larger than the 
quasiparticle coherence scale. 

Compared to the atomic limit, this approximation has a bet- 
ter behavior in the low temperature limit while still retaining a 
simplicity that allows for a completely analytic solution. The 
disadvantage of the atomic limit is that it models the zero 
temperature spectral function as an unrealistic pair of delta 
functions located at and U. The atomic limit is indeed the 
first term of the expansion in f3 J and as such is valid at high 
temperature. Our approximation instead takes into account 
the fact that the Hubbard bands are broadened by the kinetic 
term. The result is in better agreement with the DMFT data 
at temperatures lower than the limit of validity of the atomic 
limit. Our approximation eventually breaks down due to the 
fact that at low temperature the details of the shape of the Hub- 
bard bands become more relevant and the box-like model of 
the bands shows its limits. In this sense, our approximation 
should only be regarded as acceptable at high temperature. 

This form of the spectral function leads to the following ex- 
pression of the dependence of the particle number on chemical 
potential: 



n(ju, T) 



dojA(uj)f(uj) (A2) 



2k B T 




cosh 



( n-U-D \ 
\ 2k B T J 



(A3) 



In the zero-temperature limit, this expression reduces to: 

n(fi) = < —D) 

= r in + D) for — D < a < — — — 

U-A + 2D K ^ ' ' 2 

, . U — A U + A 
= 1 for < fi < 



2 . U + D. 

= 1+ 2D + U-A^-^Z— ] 

U + A n rr 

for — - — < fi < D + U 

= 2 for fi > D + U 



(A4) 



(A5) 



i.e. a function with steps at n = 0, n = 1 and n = 2 (zero 
compressibility), while n(fi) in between the steps (liquid re- 
gions) is approximated by a linear dependence on /1 (constant 
compressibility approximation). 

We regard the parameter A as a fitting parameter, function 
of coupling and temperature, and perform a least square min- 



imization of Eq. A3 to the DMFT data to determine its value. 
The result of this procedure is displayed in Fig. 26 The fit 



16 



-D 

A 



U/2-A/2 



U+D <Wf 



FIG. 25: Modelization of the spectral function in the metallic (above) 
and insulating (below) phases. In the metallic phase A = 0, while in 
the insulating phase A > 0. 
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FIG. 26: The gap parameter A(/i, T) as a function of T for different 
interaction strength U. 



as a formal limit in order to make an analytical calculation 
possible, since the approximations made in this paper are no 




3 

U/6J 



FIG. 27: Comparison of the state diagram obtained by DMFT 
(circles) and the simple analytic approximation (solid lines), [cf. 

Ref.|23") 



longer valid physically in this limit). Let us for example fo- 
cus on the crossover lines which delimit the 'Mc' region in 



Fig. 15 i within which the central occupancy is n = 1. In view 
of our approximate form of n(fi), the lower boundary of this 
regime will correspond to == (U — A)/2 and the upper 
boundary to /x = (U + A)/2. Inserting these values into 
( A7 1 and performing the integrations at T = yields the fol- 
lowing expressions for the lower and upper boundaries of the 
Mc region: 



PMc 



2V2tt U 
15 ^ 



D 



■ 2?' 2 



(A8) 



looses in quality as we reduce the temperature and the value 
of U. Restricting to the high temperature and large U region, 
the dependence of A on U and T can be separated to a good 
approximation and A(U,T) can be fit remarkably well by the 
law: 

A(U,T) ~ A + au U+a T k B T = -1.7+0.96E7 + 2Mk B T 

(A6) 

where the parameter with the largest deviation is «t which 
shows a slight increase as a function of U. 

Using the approximate expression (A3 i of n(fj,, T) for the 
homogeneous system, it is straightforward to obtain the state 
diagram of the trapped system by using the relation between 
the scaled density p and the chemical potential /iq at the center 
of the trap, which reads, in the LDA approximation: 



N 



3/2 



2vr 



^M\/Mo - V n (p,T). (A7) 



The calculations can be performed analytically in the T = 
limit (we warn the reader that T = is considered here only 
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(U- A)/D + 2 



(A9) 



Using the above determination of A, this provides an explicit 
form of the boundaries. Analytical expressions can be simi- 
larly obtained for all crossover lines in the state diagram. In 
Fig. [27] we compare these approximate analytical expressions 
to the actual lines obtained for a DMFT calculation at a low 
enough temperature T/QJ = 0.1 and find very satisfactory 
agreement. 

Finally, we note that thermodynamic quantities such as the 
double occupancy and the entropy can in principle be recon- 
structed from a given approximate expression for T, U) 
thanks to the thermodynamic (Maxwell) relations relating 
their derivatives: 



and 
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